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Objectives 

The  objective  of  our  AFOSR-sponsored  work  under  this  grant  has  been  to  advance  the  tools  that  are  required 
to  make  high-fidelity  multidisciplinary  design  of  aerospace  systems  feasible.  During  this  grant,  our  approach 
was  based  on  the  use  of  gradient-based  optimization  coupled  with  system  sensitivity  information  that  is 
inexpensively  obtained  via  adjoint  methods.  Our  effort  was  divided  into  two  main  components:  the  primary 
effort  was  focused  on  the  application  of  adjoint  techniques  to  the  problem  of  aerodynamic  shape  optimization. 
The  second  effort  developed  a  framework  to  extend  the  adjoint  approach  to  treat  additional  disciplines  in 
the  design  problem  (structures,  stability  &;  control,  propulsion,  etc.),  although  in  this  work  it  has  focused 
on  the  problem  of  aero-structural  design. 


Status  of  Effort 

The  major  objectives  stated  in  our  proposal  have  been  completed  and  achieved. 

•  During  the  first  year,  with  the  help  of  S.  Nadarajah,  a  number  of  different  studies  on  the  suitability 
of  the  continuous  and  discrete  adjoint  formulations  for  design  using  the  viscous  Reynolds-Averaged 
Navier-Stokes  equations  were  carried  out.  Our  conclusion  is  that  the  implementation  of  the  continuous 
formulation  is  more  straightforward  and  leaves  more  freedom  on  the  choices  of  discretization  of  the 
adjoint  equations.  We  found,  however,  that  the  discrete  adjoint  approach  can  be  very  helpful  in 
resolving  some  of  the  mathematical  ambiguities  that  are  left  in  the  continuous  approach  when  creating 
imposing  boundary  conditions. 

•  We  have  completed  a  new  formulation  for  the  adjoint  gradient  formulae  which  only  depends  on  the 
evaluation  of  a  surface  integral.  This  reduced  gradient  approach  has  been  tested  within  the  context  of 
structured  meshes,  but  promised  to  resolve  one  of  the  major  issues  in  the  implementation  of  adjoint 
methods  on  unstructured  meshes.  We  hope  to  pursue  this  matter  in  future  work. 

•  During  all  of  our  optimizations,  we  have  been  investigating  the  possibility  of  using  smoothing  and/or 
preconditioning  of  the  gradients  to  improve  the  convergence  and  robustness  of  the  search  procedure. 
Various  examples  are  presented  in  this  report. 

•  We  have  started  performing  planform  optimization  studies  using  the  adjoint  method.  In  the  past,  all 
of  our  work  had  fixed  the  planform  of  the  aircraft  and  had  focused  on  changes  to  the  airfoil  sections 
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along  the  span  (in  order  not  to  influence  the  other  disciplines  that  may  have  dictated  the  exact  shape 
of  the  planform.  We  are  now  in  a  position  to  obtain  accurate  gradients  of  the  cost  functions  of  interest 
to  changes  in  the  planform  shape  (sweep,  area,  taper,  etc.)  We  intend  to  continue  this  work  and  to 
include  structural  modeling  capabilities  to  ensure  that  aerodynamic  changes  in  the  planform  are  not 
detrimental  to  the  overall  performance  of  the  vehicle. 

•  We  have  developed  a  theoretical  framework  for  the  optimization  of  multi-disciplinary  systems  that 
allows  the  coupling  of  multiple  adjoint  solutions  to  provide  overall  system  sensitivities.  This  method 
is  rather  promising  in  its  application  to  high-fidelity  design. 

•  We  have  also  demonstrated  this  theoretical  framework  within  the  context  of  detailed  aero-structural 
design,  where  both  the  shape  of  the  aircraft  and  the  shape  and  material  distribution  of  the  underlying 
structure  are  simultaneously  changed  to  obtain  improvements  in  performance  measures  that  combine 
both  aerodynamic  and  structural  observations. 

In  our  opinion,  the  result  of  this  work  has  significantly  advanced  our  knowledge  in  the  application  of  adjoint 
techniques  to  both  aerodynamic  shape  optimization  and  multi-disciplinary,  aero-structural  design  problems. 
The  following  sections  summarize  the  key  results  of  this  work.  Detailed  references  are  provided  that  include 
many  more  details  on  the  derivations  and  applications  of  our  methods. 


Accomplishments 

This  section  mirrors  the  actual  work  described  above  (also  in  our  Statement  of  Objectives).  Each  section 
focuses  on  one  of  the  elements  of  our  mathematical  formulations  and  design  work. 

Continuous  vs.  Discrete  Adjoint  Method  Studies 

The  objective  of  this  study  was  to  compare  the  continuous  and  discrete  adjoint-based  automatic  aerody¬ 
namic  optimization  approaches.  Our  intention  was  to  study  the  trade-offs  between  the  complexity  of  the 
discretization  of  the  adjoint  equation  for  both  the  continuous  and  discrete  approaches,  the  accuracy  of  the 
resulting  estimates  of  the  gradients,  and  their  impact  on  the  computational  cost  to  approach  an  optimum 
solution.  For  details  on  the  derivation  and  implementation  of  the  continuous  adjoint  method,  please  refer 
to  the  cited  bibliography.  We  will  briefly  present  the  complete  formulation  of  the  discrete  adjoint  equation 
and  will  then  investigate  their  differences.  The  similarities  between  the  continuous  ’and  discrete  boundary 
conditions  are  also  explored.  The  results  demonstrate  two-dimensional  inverse  pressure  design  as  well  as  the 
accuracy  of  the  sensitivity  derivatives  obtained  from  continuous  and  discrete  adjoint-based  equations  when 
compared  to  finite-difference  gradients.  More  details  can  be  found  in  the  publications  in  the  bibliography. 

The  discrete  adjoint  equation  is  obtained  by  applying  the  control  theory  directly  to  the  set  of  discrete 
field  equations.  The  resulting  equation  depends  on  the  type  of  scheme  used  to  solve  the  flow  equations.  The 
following  discussion  uses  a  cell  centered  multigrid  scheme  with  upwind  biased  blended  first  and  third  order 
fluxes  as  the  artificial  dissipation  scheme.  A  full  discretization  of  the  equation  would  involve  discretizing 
every  term  that  is  a  function  of  the  state  vector 

nxny 

61  =  6IC  +  X)  £  (R  Hij  +  D  j)  .  (!) 
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where  6IC  is  the  discrete  cost  function,  R(w)  is  the  field  equation,  and  D(w)  is  the  artificial  dissipation 
term. 

Terms  multiplied  by  the  variation  6wij  of  the  discrete  flow  variables  are  collected  and  the  following  is 
the  resulting  discrete  adjoint  equation, 
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where 


""  ^*»i)  Ci+Jj^+2J 

+3^+I,j  W’i+i.i  “  V'l.j)  +  «i_  fjV’i-l.j 


(3) 


is  the  discrete  adjoint  artificial  dissipation  term  and  V  is  the  cell  area.  The  dissipation  coefficients  e2  and  e4 
are  functions  of  the  flow  variables,  but  to  reduce  complexity  they  are  treated  as  constants. 

In  the  case  of  an  inverse  design,  SIC  is  the  discrete  form  of  the  square  of  the  pressure  integral  around  the 
airfoil  surface.  In  contrast  to  the  continuous  adjoint,  where  the  boundary  condition  appears  as  an  update  to 
the  costate  variables  in  the  cell  below  the  wall,  the  discrete  boundary  condition  appears  as  a  source  term  in 
the  adjoint  fluxes.  At  cell  i,  2  the  adjoint  equation  is  as  follows, 


y<tk±A  — 

dt 

\  (V’i.a  -  'Pi- 1.2)  -  AT+$,2  (^i+1,2  -  'Pi, 2)] 

+  i  3-^,2)]+$, 

(4) 


where  V  is  the  cell  area,  4?  is  the  source  term  for  inverse  design, 

$  =  (—A y(ip2i,2  +  ~  (p-Pr)Asi)  Spi: 2, 


and 
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All  the  terms  in  equation  (4)  except  for  the  source  term  are  scaled  by  the  square  of  Ax.  Therefore,  as  the 
mesh  width  is  reduced,  the  terms  within  parenthesis  in  the  source  term  divided  by  A  Si  must  approach  zero  as 
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Grid  Size 

Cont. 

Disc. 

Cont-Disc 

2.397e  -  3 

9.585e  -  4 

1.724e  -  3 

2.130e  —  4 

256  x  64 

1.424e  —  3 

1.419e  -  3 

4.749e  -  5 

Table  1:  L2  norm  of  the  Difference  Between  Adjoint  and  Finite  Difference  Gradient 

the  solution  reaches  a  steady  state.  One  then  recovers  the  continuous  adjoint  boundary  condition  described 
in  earlier  works. 

If  a  first  order  artificial  dissipation  equation  is  used,  then  equation  (3)  would  reduce  to  the  term  associated 
with  e2.  In  such  a  case,  the  discrete  adjoint  equations  are  completely  independent  of  the  costate  variables  in 
the  cells  below  the  wall.  However,  if  we  use  the  blended  first  and  third  order  dissipation,  then  these  values 
are  required.  In  practice,  a  simple  zeroth  order  extrapolation  across  the  wall  produces  good  results. 

Replacing  the  inverse  design  boundary  condition  in  equation  (4)  by  the  discrete  form  of  the  cost  function 
results  in  a  discrete  adjoint  equation  for  drag  minimization. 

As  an  example  of  the  results  in  the  study,  an  inverse  design  test  case  is  briefly  shown  here.  The  target 
pressure  is  first  obtained  using  the  FL083  flow  solver  for  the  NACA  64A410  airfoil  at  a  flight  condition  of 
M  =  0.74  and  a  lift  coefficient  of  Ci  =  0.63  on  a  192  x  32  C-grid.  At  such  a  condition  the  NACA  64A410 
produces  a  strong  shock  on  the  upper  surface  of  the  airfoil,  thus  making  it  an  ideal  test  case  for  the  adjoint 
versus  finite  difference  comparison. 

The  gradient  for  the  continuous  and  discrete  adjoint  is  obtained  by  perturbing  each  point  on  the  airfoil. 
We  apply  an  implicit  smoothing  technique  to  the  gradient,  before  it  is  used  to  obtain  a  direction  of  descent 
for  each  point  on  the  surface  of  the  airfoil. 

Figures  (1),  (2),  and  (3)  exhibit  the  values  of  the  gradients  obtained  from  the  adjoint  methods  and 
finite  difference  for  various  grid  sizes.  The  circles  denote  values  that  we  obtain  by  using  the  finite  difference 
method.  The  square  represents  the  discrete  adjoint  gradient.  The  asterisk  represents  the  continuous  adjoint 
gradient.  The  gradient  is  obtained  with  respect  to  variations  in  Hicks-Henne  sine  “bump”  functions  placed 
along  the  upper  and  lower  surface  of  the  airfoil.  The  figures  only  illustrate  the  values  obtained  from  the 
upper  surface  starting  from  the  leading  edge  on  the  left  and  ending  at  the  trailing  edge  on  the  right. 

Figure  (4)  presents  the  effect  of  the  partial  discretization  of  the  flow  solver  to  obtain  the  discrete  ad¬ 
joint  equation.  Here  we  obtain  the  finite  difference  gradients  in  the  figure  without  freezing  the  dissipative 
coefficients.  A  small  discrepancy  exists  in  regions  closer  to  the  leading  edge  and  around  the  shock. 

Table  1  contains  values  of  the  L2  norm  of  the  difference  between  the  adjoint  and  finite  difference  gradients. 
The  table  illustrates  three  important  facts:  the  difference  between  the  continuous  adjoint  and  finite  difference 
gradient  is  slightly  greater  than  that  between  the  discrete  adjoint  and  finite  difference  gradient;  the  norm 
decreases  as  the  mesh  size  is  increased;  and  the  difference  between  continuous  and  discrete  adjoint  gradients 
decreases  as  the  mesh  size  is  reduced.  The  second  column  depicts  the  difference  between  the  continuous 
adjoint  and  finite  difference  gradient.  The  third  column  depicts  the  difference  between  the  discrete  adjoint 
and  finite  difference  gradients.  The  last  column  depicts  the  difference  between  the  discrete  adjoint  and 
continuous  adjoint.  As  the  mesh  size  increases  the  norms  decrease  as  expected.  Since  we  derive  the  discrete 
adjoint  by  taking  a  variation  of  the  discrete  flow  equations,  we  expect  it  to  be  consistent  with  the  finite 
difference  gradients  and  thus  to  be  closer  than  the  continuous  adjoint  to  the  finite  difference  gradient.  This 
is  confirmed  by  numerical  results,  but  the  difference  is  very  small.  As  the  mesh  size  increases,  the  difference 
between  the  continuous  and  discrete  gradients  should  decrease,  and  this  is  reflected  in  the  last  column  of 
table  1. 

In  conclusion, 

1.  The  continuous  boundary  condition  appears  as  an  update  to  the  costate  values  below  the  wall  for  a 
cell-centered  scheme,  and  the  discrete  boundary  condition  appears  as  a  source  term  in  the  cell  above 
the  wall.  As  the  mesh  width  is  reduced,  one  recovers  the  continuous  adjoint  boundary  condition  from 
the  discrete  adjoint  boundary  condition. 

2.  Discrete  adjoint  gradients  have  better  agreement  than  continuous  adjoint  gradients  with  finite  difference 
gradients  as  expected,  but  the  difference  is  generally  small. 
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Figure  1:  Adjoint  Versus  Finite  Difference 
Gradients  for  Inverse  Design  of  Korn  to 
NACA  64A410  at  Fixed  Ch 

Coarse  Grid  -  96  x  16,  M  =  0.74, 

Ci  =  0.63 


Design  Point 

Figure  3:  Adjoint  Versus  Finite  Difference 
Gradients  for  Inverse  Design  of  Korn  to 
NACA  64A410  at  Fixed  Ct. 

Fine  Grid  -  256  x  64,  M  =  0.74, 

Ci  =  0.63 


Design  Point 


Figure  2:  Adjoint  Versus  Finite  Difference 
Gradients  for  Inverse  Design  of  Korn  to 
NACA  64A410  at  Fixed  Ct. 

Medium  Grid  -  192  x  32,  M  =  0.74, 

Ci  =  0.63 


Design  Point 

Figure  4:  Adjoint  Versus  Finite  Difference 
Gradients  for  Inverse  Design  of  Korn  to 
NACA  64A410  at  Fixed  Cj.  Dissipative 
Coefficients  Not  Frozen 

Medium  Grid  -  192  x  32 ,  M  —  0.74, 

Ci  =  0.63 
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3.  As  the  mesh  size  increases,  both  the  continuous  adjoint  gradient  and  the  discrete  adjoint  gradient 
approach  the  finite  difference  gradient. 

4.  The  difference  between  the  continuous  and  discrete  gradient  reduces  as  the  mesh  size  increases. 

5.  The  cost  of  deriving  the  discrete  adjoint  is  greater. 

6.  With  our  search  procedure  as  outlined,  the  overall  convergence  of  the  objective  function  is  not  signif¬ 
icantly  affected  when  the  discrete  adjoint  gradient  is  used  instead  of  the  continuous  adjoint  gradient. 
Consequently,  we  find  no  particular  benefit  in  using  the  discrete  adjoint  method,  which  requires  greater 
computational  cost.  However,  we  believe  it  beneficial  to  use  the  discrete  adjoint  equation  as  a  guide 
for  the  discretization  of  the  continuous  adjoint  equation. 


Studies  of  Field  Integral  Elimination  in  Gradient  Computations 

A  new  continuous  adjoint  method  for  Aerodynamic  Shape  Optimization  (ASO)  using  the  Euler  and  Navier- 
Stokes  equations  has  been  developed.  Using  this  new  method  large  savings  in  the  computation  of  the 
gradients  of  aerodynamic  cost  functions  can  be  achieved  by  reducing  the  volume  integral  part  of  the  adjoint 
gradient  formula  to  a  surface  integral  over  the  configuration  of  interest.  The  new  method  offers  significant 
savings  for  three-dimensional  ASO  problems  on  general  unstructured  and  overset  meshes.  In  order  to  validate 
the  concept,  the  new  adjoint  gradient  equations  were  implemented  and  tested  for  drag  minimization  problems 
of  a  single-element  airfoil  and  a  three-dimensional  wing-fuselage  configuration.  The  accuracy  of  the  resulting 
derivative  information  for  the  two-dimensional  problem  was  investigated  by  direct  comparison  with  finite- 
difference  gradients  and  the  original  adjoint  gradients  which  include  the  volume  integral. 


The  Reduced  Gradient  Formulation 

Using  the  adjoint  formulation  that  has  been  used  in  our  previous  work,  consider  the  case  of  a  mesh  variation 
with  a  fixed  boundary.  Then, 

61  =  0, 

where  I  is  the  cost  function  of  interest.  However,  there  is  a  variation  in  the  transformed  flux, 


SFi  =  CiSw  +  SSijfj. 


Here  the  true  solution  is  unchanged, 
configuration.  Therefore 


so  the  variation  Sw  is  due  to  the  mesh  movement 


Sw  =  Vu>  ■  Sx  —  — — Sxj  (=  Sw*) 

OXj 


Sx  at  fixed  boundary 


and  since 

5-Jfl-o. 

oii 

it  follows  that 

<t>T4r{6Sijfj)dV  =  -  [  <f>TCi6w*dV.  (5) 

Jv  J'D 

A  similar  relationship  has  been  derived  in  the  general  case  with  boundary  movement  and  the  complete 
derivation  will  be  presented  in  an  upcoming  conference  paper.  Now 


<j>TSRdV 


4-Ci  (Sw  -  Sw*)  dV 
d£i 


[  4>TCi  (Sw  —  Sw*)  dB 

Jb 

f  ^-Ci(Sw-Sw*)dV. 

JV 


(6) 


Here  on  wall  boundary 


C2SW  =  SF2  -  SS2jfj- 


(7) 
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Comparison  of  Adjoint  and  Finite  Difference  Gradients 
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Figure  5:  Euler  Drag  Minimization  for  RAE2822:  Comparison  of  Original  Adjoint,  Reduced  Adjoint  and 
Finite-Difference  Gradients  Using  3  Mesh-Point  Bump  as  Design  Variable. 


Thus,  by  choosing  <j>  satisfying  the  adjoint  equation  and  the  adjoint  boundary  condition,  we  have  finally  the 
reduced  gradient  formulation  that 

SI  =  +  f  'if1'  ( 6S2jfj  +  C2$w*)  d^id^s 

J  By\> 


Gradient  comparison  for  two-dimensional  design 

Computations  were  performed  for  an  RAE  2822  airfoil  at  a  fixed  coefficient  of  lift,  Ci  —  0.6,  and  M^  =  0.75, 
using  a  C-mesh  of  size  192  x  32  for  both  the  flow  and  adjoint  solutions.  Every  mesh  point  on  the  airfoil  was 
used  as  a  design  variable.  In  this  case,  a  comparison  with  finite-difference  gradients  was  not  made,  since 
accurate  finite  difference  gradients  cannot  be  calculated  with  the  shape  discontinuity  caused  by  movement 
of  a  single  mesh  point.  The  gradients  obtained  using  both  the  original  and  new  adjoint  formulas  agree  well 
in  general.  We  believe  that  the  new  formulas  may  be  sensitive  to  details  of  their  discretization  in  regions  of 
high  curvature.  A  comparison  with  finite- difference  gradients  was  also  tried  using  local  bump  functions.  In 
order  to  avoid  geometric  discontinuities,  bumps  were  generated  such  that,  while  the  same  movement  of  each 
mesh  point  in  turn  was  made  as  before,  one- fourth  of  the  movement  was  distributed  to  its  two  neighboring 
points,  corresponding  to  local  B-splines.  Figure  5  shows  a  comparison  between  the  gradients  obtained  using 
the  three  techniques.  All  three  gradients  agree  quite  well,  although  there  is  still  a  discrepancy  near  the 
leading  edge. 


Planform  Design  Studies 

Cost  Function  for  Planform  Design 

In  order  to  design  a  high  performance  transonic  wing,  which  will  lead  to  a  desired  pressure  distribution,  and 
still  maintain  a  realistic  shape,  the  natural  choice  is  to  set 

/  —  aiC*£)  H- 02“  [  {p  —  VdfdS  +  OisCw  (®) 

*  Jb 


with 


Gw  = 


Vowing 

Qoo^ref 


(10) 
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where 

Cd  =  drag  coefficient, 

Cw  =  normalized  wing  structure  weight, 

p  =  current  surface  pressure, 

pd  =  desired  pressure, 

=  dynamic  pressure, 

Sref  =  reference  area, 

Wwing  =  wing  structure  weight,  and 

ai,  «2,  as  =  weighting  constants. 

A  practical  way  to  estimate  Wwing  is  to  use  the  so-called  Statistical  Group  Weights  Method,  which 
applies  statistical  equations  based  on  sophisticated  regression  analysis.  For  a  cargo/transport  wing  weight, 
one  can  use 

Wweight  =  0.0051  (WdgNz)°-557S™49A0'5 

+  A)01co5(A)-10ScU  (11) 


where 

A 

= 

aspect  ratio, 

Nz 

= 

ultimate  load  factor;  =  1.5  x  limit  load  factor 

Scsw 

= 

control  surface  area  (wing-mounted), 

Syj 

= 

trapezoidal  wing  area, 

tic 

= 

thickness  to  chord  ratio, 

Wdg 

= 

flight  design  gross  weight, 

A 

= 

wing  sweep,  and 

A 

= 

taper  ratio  at  25  %  MAC. 

In  addition,  if  the  wing  of  interest  is  modeled  by  five  planform  variables  such  as  root  chord  (ci),  mid-span 
chord  (02),  tip  chord  (C3),  span  (6),  and  sweepback(A),  as  shown  in  Figure  6,  it  can  be  seen  from  the  wing 
weight  formula  (11)  that  the  weight  is  estimated  to  vary  inversely  with  cos  (A),  where  A  is  the  wing  sweep. 


C3 


Figure  6:  Modeled  wing  governed  by  five  planform  variables;  root  chord  (ci),  mid-span  chord  (02),  tip  chord 
(C3),  span  (6),  and  sweepback(A). 
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Here  if  the  sweepback  is  allowed  to  vary  and  a3  is  chosen  to  be  sufficiently  large  in  the  cost  function  (9), 
we  should  expect  the  optimization  to  decrease  the  sweepback  angle  at  the  cost  of  an  increase  in  shock  drag. 

A  change  of  span  affects  the  wing  weight  of  the  wing  weight  formula  (11)  through  changes  of  the  trape¬ 
zoidal  wing  area  (Sw)j  the  aspect  ratio  (A), and  the  wing-mounted  area  (Scsw)-  From  the  wing  weight 
formula  (11),  an  increase  of  span  will  cause  SW1  A,  and  Scsw  to  increase,  resulting  in  an  increase  of  wing 
weight.  Since  induced  drag  varies  inversely  with  the  square  of  span,  if  the  span  is  allowed  to  vary  and  a3  is 
chosen  to  be  sufficiently  large,  it  is  expected  the  optimizer  to  reduce  the  span  at  the  cost  of  an  increase  in 
drag. 

Variations  of  ci,  c2,  and  c3  affect  the  wing  weight  of  the  wing  weight  formula  (11)  via  variations  of  aspect 
ratio  (A),  wing-mounted  wing  area  ( Scsw )>  trapezoidal  wing  area  (S^),  and  thickness  to  chord  ratio  ( t/c ). 
The  effect  of  an  individual  chord  change  on  the  wing  weight  is  plotted  in  Figure  7.  Figure  7  shows  that 
the  wing  weight  of  the  Statistical  Group  Weights  Method  is  a  linear  function  of  chord  length.  And  because 
the  slope  varies  along  the  span  location,  the  change  of  chord  at  different  span  location  will  affect  the  wing 
weight  differently.  With  the  same  change  in  chord  length,  the  decrease  of  the  mid-span  chord  (c2)  tends  to 
give  more  weight  reduction  than  the  others. 


Effect  of  the  planform  variables  on  the  Statistical  Group  Weights  Model 


Figure  7:  Effect  of  sweepback(A),  span  (6),  root  chord(ci),  mid-span  chord (c2),  and  tip  chord(c3)  on  the 
Statistical  Group  Weights  Method 

If  ci,  c2,  or  c3  is  allowed  to  vary,  and  again  a3  is  chosen  to  be  sufficiently  large  in  the  cost  function  (9), 
the  optimizer  would  be  expected  to  reduce  the  chord  length  with  a  consequent  increase  in  shock  drag. 

In  these  ways  the  inclusion  of  a  weight  estimate  in  the  cost  function  should  prevent  the  optimization 
from  leading  to  an  unrealistic  wing  planform,  and  yield  to  a  good  overall  performance. 
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Aerodynamic  Gradient  Calculation  for  Planform  Variables 

Gradient  information  can  be  computed  using  a  variety  of  approaches  such  as  the  finite-difference  method, 
the  complex  step  method,  and  the  automatic  differentiation.  Unfortunately,  their  computational  cost  is 
still  proportional  to  the  number  of  design  variables  in  the  problem.  In  an  optimum  transonic  wing  design, 
suppose  one  chooses  mesh  points  on  a  wing  surface  as  the  design  variables,  which  is  on  the  order  of  1000 
or  more;  it  is  impractical  to  calculate  the  gradient  using  the  methods  mentioned  earlier.  In  our  planform 
optimization,  the  design  variables  are  points  on  the  wing  surface  plus  the  planform  variables.  To  evaluate 
the  aerodynamic  gradient  with  respect  to  the  planform  variables,  since  the  number  of  planform  variables 
(five  in  this  study)  is  far  less  than  that  of  the  surface  optimization,  one  could  calculate  the  gradient  by  the 
finite-difference  method,  the  complex  step  method  or  the  automatic  differentiation.  However,  the  cost  for  the 
gradient  calculation  will  be  five  times  higher.  A  more  efficient  approach  is  to  follow  the  adjoint  formulation. 

Consider  the  aerodynamic  contribution  of  the  cost  function  (9) 


SM  dBt-  + 


v 


This  can  be  split  as 


81  =  [Iw)j  8w  +  8Iji 


8M  =  [Mw)j8w  +  8Mjr 

where  the  subscripts  I  and  II  are  used  to  distinguish  between  the  contributions  associated  with  variation  of 
the  flow  solution  8w  and  those  associated  with  the  metric  variations  <55.  Thus  [Mw]i  represents  with 
the  metrics  fixed.  Note  that  SR  is  intentionally  kept  unsplit  for  programming  purposes.  If  one  chooses  ip  as 
ip *  where  ip *  satisfies 


6I(w,  S) 


6I(S) 

J  SMndB(  +  J  ip* T6RdV t 

B  V 

J^SMnAB  +  ^r^R 

B  V 

£ SMnAB  +  (R\s+ss  -  R\s) , 

B  V 


where  R\$  and  R\s+ss  are  volume  weighted  residuals  calculated  at  the  original  mesh  and  at  the  mesh 
perturbed  in  the  design  direction. 

Provided  that  ip*  has  already  been  calculated  and  R  can  be  easily  calculated,  the  gradient  of  the  planform 
variables  can  be  computed  effectively  by  first  perturbing  all  the  mesh  points  along  the  direction  of  interest. 
For  example,  to  calculate  the  gradient  with  respect  to  the  sweepback,  move  all  the  points  on  the  wing  surface 
as  if  the  wing  were  pushed  backward  and  also  move  all  other  associated  points  in  the  computational  domain 
to  match  the  new  location  of  points  on  the  wing.  Then  re-calculate  the  residual  value  and  subtract  the 
previous  residual  value  from  the  new  value  to  form  A R.  Finally,  to  calculate  the  planform  gradient,  multiply 
A R  by  the  costate  vector  and  add  the  contribution  from  the  boundary  terms. 

This  way  of  calculating  the  planform  gradient  exploits  full  benefit  of  knowing  the  value  of  adjoint  variables 
ip*  with  no  extra  cost  of  flow  or  adjoint  calculations. 


Choice  of  Weighting  Constants 

The  choice  of  aq  and  o3  greatly  affects  the  optimum  shape.  An  intuitive  choice  of  and  a3  can  be  made  by 
considering  the  problem  of  maximizing  range  of  an  aircraft.  The  simplified  range  equation  can  be  expressed 
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R=V-Llog^l 

CD  yW2 


where 

C  =  Specific  Fuel  Consumption, 

D  =  Drag, 

L  =  Lift, 

R  =  Range, 

V  =  Aircraft  velocity, 

W\  —  Take  off  weight,  and 

W2  =  Landing  weight. 

If  one  takes 


Wi  =  We  +  Wf  =  fixed 
W2  =  We 


W*  =  Gross  weiSht  of  the  airplane  without  fuel,  then  ^  variation  of  the  weight  can  be  expressed  as 
Wf  =  Fuel  weight, 

SW2  =  SWe. 

With  fixed  ^  ,  Wu  and  L,  the  variation  of  R  can  be  stated  as 

« -  £('(£)<♦£«)) 


Vf.(L\,  W,L(  W, 
c(‘{d)  109  W2  +  DS{,0aW2 
V  f  6DL  Wi  L6W2\ 
C\  D  Dl°9W2  DW2) 


V  L,  Wi  ( SD  1  6W2 

CD109W2{d  +  log&  w2 


If  we  minimize  the  cost  function  defined  as 


SCD  1  SW2 
Cd  +  log IV2 


5Cd  1  SCw2 

Cd  log ^  Cw2 


I  =  Cd  +  dCwy 


where  a  is  the  weighting  multiplication,  then  choosing 


cw2log ^ 


corresponds  to  maximizing  the  range  of  the  aircraft. 


Design  Cycle 

The  design  cycle  starts  by  first  solving  the  flow  field  until  at  least  a  4  order  of  magnitude  drop  in  the 
residual.  The  flow  solution  is  then  passed  to  the  adjoint  solver.  Second,  the  adjoint  solver  is  run  to  calculate 
the  costate  vector.  Iteration  continues  until  at  least  a  4  order  of  magnitude  drop  in  the  residual.  The  costate 
vector  is  passed  to  the  gradient  module  to  evaluate  the  aerodynamic  gradient.  Then,  the  structural  gradient 
is  calculated  and  added  to  the  aerodynamic  gradient  to  form  the  overall  gradient.  The  steepest  descent 
method  is  used  with  a  small  step  size  to  guarantee  that  the  solution  will  converge  to  the  optimum  point. 
The  design  cycle  is  shown  in  Figure  8. 
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Figure  8:  Design  cycle 


Applications  to  Aero- Structural  Optimization 

During  the  course  of  this  grant,  we  have  developed  and  implemented  a  high-fidelity  aero-structural  opti¬ 
mization  framework  that  allows  the  optimization  of  aircraft  configurations  described  by  the  full  Euler  or 
Navier-Stokes  equations  and  a  complete  linear  finite  element  model  of  the  structure.  The  framework  has 
been  recently  completed  and  uses  a  coupled  aero-structural  adjoint  procedure  to  compute  coupled  sensitiv¬ 
ities  of  both  aerodynamic  and  structural  cost  functions  to  design  variables  that  parameterize  the  shape  of 
both  the  outer  mold  line  and  the  structure  itself.  The  accuracy  of  the  sensitivities  has  been  demonstrated 
and  their  use  in  realistic  design  has  also  been  presented.  In  order  to  reduce  the  cost  of  computation  of  the 
sensitivities  of  the  finite  element  stresses,  Kreisselmeier-Steinhauser  (K-S)  functions  are  used  to  lump  the 
constraints  in  the  problem. 

With  a  suitable  definition  of  the  design  variables  and  constraints,  we  are  seeking  to  solve  the  following 
aircraft  design  optimization  problem: 

minimize  I  —  ctCo  +  fiW 
XA,Xs€l&n 
subject  to  Cl  =  Clt 
KS  >  0 
xs>xsmin- 

with  the  use  of  the  coupled  adjoint  procedure  to  calculate  the  sensitivities  required  by  the  optimizer. 


Aero-Structural  Sensitivity  Analysis 


For  the  case  of  aero- structural  sensitivities,  we  have  coupled  aerodynamic  ( Ra )  and  structural  ( Rs )  governing 
equations,  and  two  sets  of  state  variables:  the  flow  state  vector,  w,  and  the  vector  of  structural  displacements, 
u.  In  the  following  expressions,  we  split  the  vectors  of  residuals,  states  and  adjoints  into  two  smaller  vectors 
corresponding  to  the  aerodynamic  and  structural  systems,  i.e. 


R  = 


'  Ra  ' 

w 

ill  — 

IpA 

.  Rs  . 

>  y  — 

u 

>  w  — 

ips  . 

Figure  9  shows  a  diagram  representing  the  coupling 
equation  for  an  aero-structural  system  can  be  written 


(13) 

new  notation,  the  adjoint 


as 


r  Ma 

dw 


i 

du  J 


iT 


IpA 

IpS 


aj_  -] 

du  . 


(14) 
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Figure  9:  Schematic  Representation  of  the  Aero-Structural  Governing  Equations. 

In  addition  to  the  diagonal  terms  of  the  matrix  that  appear  when  we  solve  the  single  discipline  adjoint 
equations,  we  also  have  off-diagonal  terms  expressing  the  sensitivity  of  one  discipline  to  the  state  variables 
of  the  other.  The  residual  sensitivity  matrix  in  this  equation  is  identical  to  that  of  the  Global  Sensitivity 
Equations  (GSE)  introduced  by  Sobieski. 

Since  the  factorization  of  the  full  matrix  in  the  system  of  equations  (14)  would  be  extremely  costly,  our 
approach  uses  an  iterative  solver,  much  like  the  one  used  for  the  aero-structural  solution,  where  the  adjoint 
vectors  are  lagged  and  the  two  different  sets  of  equations  are  solved  separately.  For  the  calculation  of  the 
adjoint  vector  of  one  discipline,  we  use  the  adjoint  vector  of  the  other  discipline  from  the  previous  iteration, 
i.e.,  we  solve 

9Ra 

dw 


,  01  ORs  *  7 


(16) 


where  ip  a  and  tps  are  the  lagged  aerodynamic  and  structural  adjoint  vectors.  The  final  result  given  by  this 
system,  is  the  same  as  that  of  the  original  coupled- adjoint  equations  (14).  We  call  this  the  Lagged- Coupled 
Adjoint  (LCA)  method  for  computing  sensitivities  of  coupled  systems.  Note  that  these  equations  look  like 
the  single  discipline  adjoint  equations  for  the  aerodynamic  and  the  structural  solvers,  with  the  addition  of 
forcing  terms  in  the  right-hand-side  that  contain  the  off-diagonal  terms  of  the  residual  sensitivity  matrix. 
Note  also  that,  even  for  more  than  two  disciplines,  this  iterative  solution  procedure  is  nothing  but  the 
well-known  Block- Jacobi  method. 

As  noted  previously,  dRs/du  =  K  for  a  linear  structural  solver.  Since  the  stiffness  matrix  is  symmetric 
(Kt  =  K)  the  structural  equations  (16)  are  self-adjoint.  Therefore,  the  structural  solver  can  be  used  to 
solve  for  the  structural  adjoint  vector,  ^5,  by  using  the  pseudo-load  vector  given  by  the  right-hand-side  of 
equation  (16). 

Once  both  adjoint  vectors  have  converged,  we  can  compute  the  final  sensitivities  of  the  objective  function 


by  using 


dl  dl  ,t9Ra  .  ,t9Rs 


(17) 


which  is  the  coupled  version  of  the  total  sensitivity  equation.  In  order  to  solve  the  aircraft  optimization 
problem  we  proposed  earlier  on,  we  also  need  sensitivities  of  the  structural  weight  with  respect  to  the 
design  variables.  Since  the  aero-structural  coupling  does  not  involve  the  weight,  these  sensitivities  are  easily 
computed. 

Figure  10  below  shows  the  coupled  aero-structural  sensitivities  of  the  Cd  with  respect  to  a  number  of 
aerodynamic  shape  variables  on  the  upper  surface  of  a  supersonic  business  jet  wing.  The  agreement  of  our 
method,  in  comparison  with  results  produced  with  the  complex-step  formulation  is  quite  good.  Similar 
agreement  is  found  for  structural  cost  functions  and  structural  design  variables. 
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Figure  10:  Sensitivities  of  the  Drag  Coefficient  with  Respect  to  Shape  Perturbations. 

Aero-Structural  Design 

The  objective  in  this  optimization  is  to  solve  the  design  problem  that  we  previously  described,  i.e., 

minimize  I  =  olCd  +  (3W 

xeRn 

subject  to  Cl  =  Clt 
KS  >  0 

®  ^  ®min* 

In  our  example  the  value  of  Cd  corresponds  to  that  of  the  cruise  condition,  which  has  a  target  lift 
coefficient  of  0.1.  The  structural  stresses,  in  the  form  of  the  KS  function,  correspond  to  a  single  maneuver 
condition,  for  which  CxT  =0.2. 

All  optimization  work  is  carried  out  using  the  nonlinear  constrained  optimizer  NPSOL,  Euler  calculations 
are  performed  on  a  wing-body  36-block  mesh  that  is  constructed  from  the  decomposition  of  a  193  x  33  x  49 
C-H  mesh.  During  the  process  of  optimization,  all  flow  evaluations  are  converged  to  5.3  orders  of  magnitude 
of  the  average  density  residual  and  the  Cl  constraint  is  satisfied  within  1CT6. 

In  order  to  parameterize  the  shape  of  the  aircraft,  we  have  chosen  sets  of  design  variables  that  apply  to 
both  the  wing  and  the  fuselage.  The  wing  shape  is  modified  by  the  design  optimization  procedure  at  six 
defining  stations  uniformly  distributed  from  the  side-of-body  to  the  tip  of  the  wing.  The  shape  modifications 
of  these  defining  stations  are  linearly  lofted  to  a  zero  value  at  the  previous  and  next  defining  stations.  On 
each  defining  station,  the  twist,  the  leading  and  trailing  edge  camber  distributions,  and  five  Hicks-Henne 
bump  functions  on  both  the  upper  and  lower  surfaces  are  allowed  to  vary.  The  leading  and  trailing  edge 
camber  modifications  are  not  applied  at  the  first  defining  station.  This  yields  a  total  of  76  OML  design 
variables  on  the  wing.  Planform  modifications,  which  are  permitted  by  our  software,  were  not  used  in  the 
present  calculations.  Planform  optimization  is  only  meaningful  if  additional  disciplines  and  constraints  are 
taken  into  account. 

The  shape  of  the  fuselage  is  parameterized  in  such  a  way  that  its  camber  is  allowed  to  vary  while  the  total 
volume  remains  constant.  This  is  accomplished  with  9  bump  functions  evenly  distributed  in  the  streamwise 
direction  starting  at  the  10%  fuselage  station.  Fuselage  nose  and  trailing  edge  camber  functions  are  added 
to  the  fuselage  camber  distribution  in  a  similar  way  to  what  was  done  with  the  wing  sections. 

The  structural  sizing  is  accomplished  with  10  design  variables,  which  correspond  to  the  skin  thicknesses 
of  the  top  and  bottom  surfaces  of  the  wing.  Each  group  is  formed  by  the  plate  elements  located  between  two 
adjacent  ribs.  All  structural  design  variables  are  constrained  to  exceed  a  specified  minimum  gauge  value. 

The  complete  configuration  is  therefore  parameterized  with  a  total  of  97  design  variables.  As  mentioned 
in  an  earlier  section,  the  cost  of  aero-structural  gradient  information  using  our  coupled-adjoint  method  is 
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Figure  11:  Convergence  history  of  the  aero-structural  optimization. 


effectively  independent  of  the  number  of  design  variables:  in  more  realistic  full  configuration  test  cases  that 
we  are  about  to  tackle,  500  or  more  design  variables  will  be  necessary  to  describe  the  shape  variations  of  the 
configuration  (including  nacelles,  diverters,  and  tail  surfaces)  and  the  sizing  of  the  structure. 

The  initial  application  of  our  design  methodology  to  the  aero-structural  design  of  a  supersonic  business  jet 
is  simply  a  proof- of-concept  problem  meant  to  validate  the  sensitivities  obtained  with  our  method.  Current 
work  is  addressing  the  use  of  multiple  realistic  load  conditions,  dynamic  loads,  aeroelastic  constraints,  and 
the  addition  of  diverters,  nacelles,  and  empennage. 

In  the  present  design  case,  we  use  a  =  104,  =  3.226  x  1CT3.  Note  that  the  scalars  that  multiply  the 

structural  weight,  W,  and  the  coefficient  of  drag,  Cd ,  reflect  the  correct  trade-off  between  drag  and  weight 
that  was  previously  mentioned,  i.e.  that  one  count  of  drag  is  worth  310  pounds  of  weight. 

Figure  11  shows  the  evolution  of  this  aero-structural  design  case  for  successive  major  design  iterations. 
The  figure  shows  the  values  of  the  coefficient  of  drag  (in  counts),  the  wing  structural  weight  (in  lbs),  and  the 
value  of  the  KS  function.  Note  that  the  structural  constraints  are  satisfied  when  the  KS  function  is  positive. 
Because  of  the  approximate  nature  of  the  KS  function,  all  structural  constraints  may  actually  be  satisfied 
for  small  but  negative  values  of  the  KS  function. 

The  baseline  design  is  feasible,  with  a  cruise  drag  coefficient  of  74.04  counts  and  a  structural  weight 
of  9,285  lbs.  The  KS  function  is  slightly  positive  indicating  that  all  stress  constraints  are  satisfied  at  the 
maneuver  condition.  In  the  first  two  design  iterations,  the  optimizer  takes  large  steps  in  the  design  space, 
resulting  in  a  drastic  reduction  in  both  Cd  and  W.  However,  this  also  results  in  a  highly  infeasible  design  that 
exhibits  maximum  stresses  that  have  a  value  of  2.1  times  the  yield  stress  of  the  material.  After  these  initial 
large  steps,  the  optimizer  manages  to  decrease  the  norm  of  the  constraint  violation.  This  is  accomplished 
by  increasing  the  structural  skin  thicknesses  while  decreasing  the  airfoil  thicknesses,  resulting  in  a  weight 
increase  and  a  further  reduction  in  drag.  Towards  major  iteration  10,  there  is  no  visible  progress  for  several 
iterations  while  the  design  remains  infeasible.  In  iteration  13,  a  large  design  step  results  in  a  sudden  increase 
in  feasibility  accompanied  by  an  equally  sudden  increase  in  Cd  ■  The  optimizer  has  established  that  the  best 
way  of  obtaining  a  feasible  design  is  to  increase  the  wing  thickness  (with  the  consequent  increases  in  Cd  and 
weight)  and  the  structural  thicknesses.  From  that  point  on,  the  optimizer  rapidly  converges  to  the  optimum. 
After  43  major  iterations,  the  KS  constraint  is  reduced  to  (9(10-4)  and  all  stress  constraints  are  satisfied. 
The  aero-structurally  optimized  result  has  CD  =  0.006922  and  a  total  wing  structure  weight  of  5,546  lbs. 

Visualizations  of  the  baseline  and  optimized  configurations  are  shown  in  Figures  12  and  13.  Measures  of 
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Figure  12:  Baseline  configuration  for  the  supersonic  business  jet  showing  surface  densities  at  the  cruise 
condition  and  structural  stresses  at  the  maneuver  condition.  The  density  is  normalized  by  the  freestream 
value  and  the  von  Mises  stresses  are  normalized  by  the  material  yield  stress. 
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performance  and  feasibility  are  presented  in  the  first  section  of  Table  2.  The  left  halves  of  Figures  12  and  13 
show  the  surface  density  distributions  with  the  corresponding  structural  deflections  at  the  cruise  condition 
for  both  the  initial  and  optimized  designs.  The  right  halves  show  exploded  views  of  the  stress  distributions 
on  the  structure  (spar  caps,  spar  shear  webs,  and  skins,  from  top  to  bottom)  at  the  Cl  =  0.2  maneuver 
condition.  From  these  Figures  one  can  appreciate  that  not  only  have  the  surface  density  distributions 
changed  substantially  at  the  cruise  point,  but  so  have  the  element  stresses  at  the  maneuver  condition.  In 
fact,  as  expected  from  a  design  case  with  a  single  load  condition,  the  optimized  structure  is  exhibits  stresses 
much  closer  to  the  yield  stress,  except  in  the  outboard  sections  of  the  wing,  where  the  minimum  gauge 
constraints  are  active.  It  is  also  worth  noting  that  about  half  of  the  improvement  in  the  Cd  of  the  optimized 
configuration  results  from  drastic  changes  in  the  fuselage  shape:  both  front  and  aft  camber  have  been  added 
to  distribute  the  lift  more  evenly  in  the  streamwise  direction  in  order  to  reduce  the  total  lift-dependent  wave 
drag. 

A  total  of  50  major  design  iterations  including  aero-structural  analyses,  coupled  adjoint  solutions,  gradient 
computations,  and  line  searches  were  performed  in  approximately  20  hours  of  wall  clock  time  using  18 
processors  of  an  SGI  Origin  3000  system  (R12000,  400  MHz  processors).  Since  these  are  not  the  fastest 
processors  currently  available  we  feel  confident  that  much  larger  models  can  be  optimized  with  overnight 
turnaround  in  the  near  future. 

Comparison  with  Sequential  Optimization 

The  usefulness  of  a  coupled  aero-structural  optimization  method  can  only  be  measured  by  comparing  with 
the  results  obtained  using  current  state-of-the-art  practices.  In  the  case  of  aero-structural  design,  the  typical 
approach  is  to  carry  out  aerodynamic  shape  optimization  with  artificial  airfoil  thickness  constraints  meant  to 
represent  the  effect  of  the  structure,  followed  by  structural  optimization  with  a  fixed  OML.  It  is  well  known 
that  sequential  optimization  cannot  be  guaranteed  to  converge  to  the  true  optimum  of  a  coupled  system.  In 
order  to  determine  the  difference  between  the  optima  achieved  by  fully-coupled  and  sequential  optimizations, 
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Figure  13:  Optimized  configuration  for  the  supersonic  business  jet. 
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we  have  also  carried  out  one  cycle  of  sequential  optimization  within  our  analysis  and  design  framework. 

To  prevent  the  optimizer  from  thinning  the  wing  to  an  unreasonable  degree  during  the  aerodynamic 
shape  optimization,  5  thickness  constraints  are  added  to  each  of  the  6  defining  stations  for  a  total  of  30 
linear  constraints.  These  constraints  are  such  that,  at  the  points  where  they  are  applied,  the  wing  box  is  not 
allowed  to  get  any  thinner  than  the  original  design. 

After  the  process  of  aerodynamic  shape  optimization  is  completed,  the  initial  Cp  has  decreased  to 
0.006992,  as  shown  in  the  lower  portion  of  Table  2.  After  fixing  the  OML,  structural  optimization  is 
performed  using  the  maneuver  loads  for  the  baseline  configuration  at  Cl  —  0.2.  The  structural  optimization 
process  reduces  the  weight  of  the  wing  structure  to  6,567  lbs. 

We  can  now  compare  the  results  of  the  fully  coupled  optimization  in  the  previous  section  and  the  outcome 
of  the  process  of  sequential  optimization.  The  differences  are  clear:  the  coupled  aero-structural  optimization 
was  able  to  achieve  a  design  with  a  range  of  7,361  nm,  which  is  224  nm  higher  than  that  obtained  from  the 
sequential  optimization. 

Finally,  note  that  since  sequential  optimization  neglects  the  aero-structural  coupling  in  the  computation 
of  maneuver  loads,  there  is  no  guarantee  that  the  resulting  design  is  feasible.  In  fact,  the  aero-structural 
analysis  shows  that  the  value  of  the  KS  function  is  slightly  negative. 
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Cd  (counts) 

KS 

U'max/G’y 

ZFW  (lbs) 

Range  (nm) 

Baseline 

73.95 

1.15  x  1CT1 

0.87 

47,500 

6,420 

Integrated  optimization 

69.22 

-2.68  x  10-4 

0.98 

43,761 

7,361 

Sequential  optimization 

Aerodynamic  optimization 

Baseline 

74.04 

Optimized 

69.92 

Structural  optimization 

Baseline 

1.02  x  10_1 

0.89 

47,500 

Optimized 

1.45  x  10~8 

0.98 

44,782 

Aero-structural  analysis 

69.92 

-9.01  x  10“3 

0.99 

7,137 

Table  2:  Comparison  between  the  integrated  and  sequential  approaches  to  aero-structural  optimization. 
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